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We study the transition from inspiral to plunge in general relativity by computing gravitational 
waveforms of non-spinning, equal-mass black-hole binaries. We consider three sequences of simula- 
tions, starting with a quasi-circular inspiral completing 1.5, 2.3 and 9.6 orbits, respectively, prior to 
coalescence of the holes. For each sequence, the binding energy of the system is kept constant and 
the orbital angular momentum is progressively reduced, producing orbits of increasing eccentricity 
and eventually a head-on collision. We analyze in detail the radiation of energy and angular mo- 
mentum in gravitational waves, the contribution of different multipolar components and the final 
spin of the remnant, comparing numerical predictions with the post-Newtonian approximation and 
with extrapolations of point-particle results. We find that the motion transitions from inspiral to 
plunge when the orbital angular momentum L = Lcrit — 0.8M^. For L < Lcrit the radiated energy 
drops very rapidly. Orbits with L ~ Lcrit produce our largest dimensionless Kerr parameter for 
the remnant, j = J/M'^ ~ 0.724 ± 0.13 (to be compared with the Kerr parameter j ~ 0.69 result- 
ing from quasi-circular inspirals). This value is in good agreement with the value of 0.72 reported 
in [l[. These conclusions are quite insensitive to the initial separation of the holes, and they can 
be understood by extrapolatingpoint particle results. Generalizing a model recently proposed by 
Buonanno, Kidder and Lehner Q] to eccentric binaries, we conjecture that (1) j ~ 0.724 is close to 
the maximal Kerr parameter that can be obtained by any merger of non-spinning holes, and (2) no 
binary merger (even if the binary members are extremal Kerr black holes with spins aligned to the 
orbital angular momentum, and the inspiral is highly eccentric) can violate the cosmic censorship 
conjecture. 

PACS numbers: 04.25.dg, 04.25.Nx, 04.30.Db, 04.70.Bw 



I. INTRODUCTION 



The research area of gravitational wave (GW) physics has reached a very exciting stage, both experimentally 
and theoretically. Earth-based laser-interferometric detectors, including LIGO GEO600 Q and TAMA are 
collecting data at design sensitivity, searching for GWs in the frequency range ~ 10 — lO'^ Hz. VIRGO [gI] should 
reach design sensitivity within one year, and the space-based interferometer LISA is expected to open an observational 
window at low frequencies (~ 10~^ — 10^^ Hz) within the next decade 0. 

The last two years have also seen a remarkable breakthrough in the simulation of the strongest expected GW 
sources, the inspiral and coalescence of black-hole binaries [8|, |9|. |ip|. Several groups have now generated independent 
numerical codes for such simulations [ll|i Eli IHi Elj Il7l. ll8l, |l9l| and studied various aspects of binary black hole 



mergers. In the context of analyzing the resulting gravitational waveforms, these include in particular the comparisons 



of numerical results with post-Newtonian (PN) predictions [20|, |21|, 
radiation 




multipolar analyses of the emitted 
ISOj l and gravitational wave emission 



|2J, |26| , the use of numerical waveforms in data analysis [27 

from systems of three black holes 3l|. 

Despite this progress, a comprehensive analysis of binary black hole inspirals remains a daunting task, mainly 
because of the large dimensionality of the parameter space. In geometrical units, the total mass of the binary is just 
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an overall scale factor. The source parameters to be explored by numerical simulations (sometimes called "intrinsic" 
parameters in the GW data analysis literature) include the mass ratio q — M2/M1, the eccentricity e of the orbit and 
six parameters for the magnitude of the individual black hole spins and their direction with respect to the binary's 
orbital angular momentum. 

In this paper we present results from numerical simulations of non-spinning, equal-mass black-hole binaries, and 
we focus on the effect of the orbital eccentricity on the merger waveforms. We consider three sequences, starting with 
quasi-circular inspirals that complete ~ 1.5, ~ 2.3 and ~ 9.6 orbits, respectively, prior to coalescence of the holes. 
By fixing the binding energy of the system and progressively reducing the orbital angular momentum, we produce 
a sequence of orbits of increasing eccentricity and eventually a head-on collision. For each of these simulations we 
analyze in detail the radiation of energy and angular momentum in GWs, the contribution of different multipolar 
components and the final spin of the remnant, comparing numerical predictions with the PN approximation and with 
extrapolations of point-particle results. 

Non-eccentric inspirals are usually considered the most interesting cases for GW detection. For an isolated binary 
evolving under the effect of gravitational radiation reaction, the eccentricity decreases by roughly a factor of 3 when 
the orbital semimajor axis is halved fs^]. For most conceivable formation mechanisms of solar-mass black hole 
binaries, the orbit will usually be circular by the time the GW signal enters the best-sensitivity bandwidth of Earth- 
based interferometers. However, we wish to stress that our simulations could be of interest for GW detection. For 
example, according to some astrophysical scenarios, eccentric binaries may be potential GW sources for Earth-based 
detectors. In globular clusters, the inner binaries of hierarchical triplets undergoing Kozai oscillations can merge 
under gravitational radiation reaction, and ^ 30% of these systems can have eccentricity ~ 0.1 when GWs enter the 
detectors' most sensitive bandwidth at ^ 10 Hz [33j. Massive black hole binaries to be observed by LISA could also 
have significant eccentricity in the last year of inspiral. Recent simulations using smoothed particle hydrodynamics 
follow the dynamics of binary black holes in massive, rotationally supported circumnuclear discs [3^. l35l. [36| . In these 
simulations, a primary black hole is placed at the center of the disc and a secondary black hole is set initially on an 
eccentric orbit in the disc plane. By using the particle splitting technique, the most recent simulations follow the 
binary's orbital decay down to distances ~ 0.1 pc. Dynamical friction is found to circularize the orbit if the binary 
corotates with the disc [35]. However, if the orbit is counterrotating with the disc the initial eccentricity does not seem 
to decrease, and black holes may still enter the GW emission phase with high eccentricity [3^ . 

Complementary studies show that eccentricity evolution may still occur in later stages of the binary's life, because 
of close encounters with single stars and/or gas-dynamical processes. Three-body encounters with background stars 
have been studied mainly in spherical backgrounds. These studies find that stellar dynamical hardening can lead to 
an increase of the eccentricity, acting against the circularization driven by the large-scale action of the gaseous and/or 
stellar disc, possibly leaving the binary with non-zero eccentricity when gravitational radiation reaction becomes 
dominant [33, [H, [s^, HO, [4l[. It has also been suggested that the gravitational interaction of a binary with a 
circumbinary gas disc could increase the binary's eccentricity. The transition between disc-driven and GW-driven 
inspiral can occur at small enough radii that a small but significant eccentricity survives, typical values being e ~ 0.02 
(with a lower limit e ~ 0.01) one year prior to merger (cf. Fig. 5 of [i^l). If the binary has an "extreme" mass ratio 
q < 0.02 the residual eccentricity predicted by this scenario can be considerably larger, e > 0.1. Numerical simulations 
should be able to test these predictions in the near future. As shown by Sopuerta, Yunes and Laguna, eccentricity 
could significantly increase the recoil velocity resulting from the merger of non-spinning black-hole binaries [43|. 

Independently of the presence of eccentricity in astrophysical binary mergers, the problem we consider here has 
considerable theoretical interest. Our simulations explore the transition between gravitational radiation from a quasi- 
circular inspiral (the expected final outcome in most astrophysical scenarios) and the radiation emitted by a head-on 
collision, where the binary has maximal symmetry. Our work should provide some guidance for analytical studies 
of the "transition from inspiral to plunge" . The first analytical study of this problem in the context of PN theory 
was carried out by Kidder, Will and Wiseman [ij. The transition between the adiabatic phase and the plunge was 
studied in [i^ using nonperturbative resummed estimates of the damping and conservative parts of the two-body 
dynamics, i.e. the so-called "Effective One Body" (EOB) model. Ori and Thorne 'i^ provided a semi-analytical 
treatment of the transition in the extreme mass ratio limit. Waveforms comprising inspiral, merger and ringdown 
for comparable-mass bodies have also been produced using the EOB model (see eg. [43| for extensions of the original 
model to spinning binaries and for references to previous work). Preliminary comparisons of EOB and numerical 
relativity waveforms showed that improved models of ringdown excitation [23, 28, 48] or additional phenomenological 
terms in the EOB effective potential [i^ are needed to achieve acceptable phase differences between the numerical 
and analytical waveforms. 

Our study is complementary to Ref. (soj . that considered sequences of eccentric, equal-mass, non-spinning binary 
black hole evolutions around the "threshold of immediate merger" : a region of parameter space separating binaries 
that quickly merge to form a final Kerr black hole from those that do not merge in a short time. Similar scenarios 
have also been studied in Ref. [5l| , with particular regard to the maximal spin of the final hole generated in this way. 
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The universality of the gravitational wave signal during the merger was analysed in Ref. where it was pointed 
out that binaries largely circularize after about 9 orbits when starting with eccentricities below about 0.4. The first 
comparison between numerical evolutions of eccentric binaries with post-Newtonian predictions was presented in [s^] . 

Our focus in this work is on the high-eccentricity region of the parameter space, which always leads to merger. In 
particular, the near-head-on limit of our study is of interest as a first step to compute the energy loss and production 
cross-section of mini-black holes in TeV-scale gravity scenarios (possibly at the upcoming LHC [E^l), and trans- 
Planckian scattering in general [s^ [55| . Present semi-analytical techniques (including a trapped surface search in 
the union of Aichelburg-Sexl shock waves, close-limit approximation calculations and perturbation theory) only give 
rough estimates of the emitted energy and production cross-section [5^ and do not provide much insight into the 
details of the process (but see [131 fo'' ^ ^i'^* numerical investigation). 

Our main finding is that, for all sequences we studied, the motion radically changes character when the black holes' 
orbital angular momentum L ^ icrit — 0.8M^, turning from an eccentric inspiral into a plunge. In particular, for 
L < Lcrit we observe that: 

• The number of orbits iV-i^avcs (as estimated using the gravitational wave cycles) or A^punc (as computed from the 
punctures' trajectories) becomes less than one, so the motion effectively turns into a plunge (see Table U and 
Fig. m below); 

• The energy emission starts decreasing exponentially (Fig. [7]); 

• PN-based eccentricity estimates yield meaningless results (Table |l|; 

• The polarization becomes linear rather than circular (Fig. ^ ; 

• The final angular momentum starts decreasing, rather than increasing, as P and L decrease (Fig. 

Binary mergers with L ~ Lcrit are those producing the largest Kerr parameter for the final black hole observed in 
our simulations, jfin — 0.724. One is led to suspect that for maximally spinning holes having spins aligned with the 
orbital angular momentum, a large orbital eccentricity may lead to violations of the cosmic censorship conjecture. 
Using arguments based on the extrapolation of point-particle results (see also pi]), we conjecture that (1) the maximal 
Kerr parameter that can be obtained by any merger of non-spinning holes is not much larger than j ~ 0.724, and (2) 
cosmic censorship will not be violated as a result of any merger, even in the presence of orbital eccentricity. Further 
numerical simulations are needed to confirm or disprove these conjectures. 

The paper is organized as follows. We begin in Sec.[n]discussing to what extent the Newtonian concept of eccentricity 
can be generalized to characterize orbiting binaries in general relativity. For this purpose, we introduce and compare 
various PN estimates of the orbital eccentricity, and we show that these eccentricity estimates break down when the 
motion turns from inspiral to plunge. Sec. IIIII is a brief introduction to the numerical code used for the simulations. 
After a discussion of the choice of initial data and of the code's accuracy, we show how reducing the orbital angular 
momentum affects the gravitational waveforms, the puncture trajectories and the polarization of the waves. In Sec. lIVI 
we study the multipolar energy distribution of the radiation and the angular momentum of the final Kerr black hole. 
In Sec.|V]we show that the salient features of our simulations can be understood using extrapolations of point-particle 
results. Sec. IVII is devoted to fits of the ringdown waveform and to estimates of the energy radiated in ringdown 
waves. We conclude by considering possible future extensions of our investigation. 

II. POST-NEWTONIAN ESTIMATES OF THE ECCENTRICITY 

In Newtonian dynamics, the shape of a binary's orbital configuration is determined by two parameters, the semi- 
major axis and the eccentricity. These parameters are intimately tied to the binding energy and orbital angular 
momentum of the binary and our construction of sequences of binaries with increasing eccentricity is based on this 
Newtonian intuition. Specifically, we fix the binding energy of the system, progressively reduce the orbital angular 
momentum and thus produce a sequence of orbits of increasing eccentricity. Before doing so, however, we need to 
address a conceptual difficulty, namely, how to quantify eccentricity in general relativity. 

It turns out, unfortunately, that there exists no unique, unambiguous definition of eccentricity in fully non-linear 
general relativity. For this reason, in the following we will use PN arguments to quantify the initial eccentricity (or 
rather, eccentricities) of the simulations. We will consider in detail two different generalizations of the Newtonian 
eccentricity: the 3PN extension [ssj of a quasi-Keplerian parametrization originally proposed by Damour and Deruelle 
[sgj , and a definition in terms of observable quantities recently introduced by Mora and Will ^] . 
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FIG. 1; The PN eccentricity parameters et, and for an equal mass binary with binding energy _Eb are shown as functions 
of the orbital angular momentum LjM^ for ADM-type coordinates and harmonic coordinates. The left panel shows the result 
for a binding energy corresponding to our sequence 1, the right panel that obtained for a much smaller binding energy. 



A. Quasi-Keplerian parametrization 

A quasi-Keplerian parametrization of eccentric orbits of objects with mass M\ and Mi has been derived at IPN 
order in harmonic coordinates by Damour and Deruelle [5§|, extended to 2PN order in ADM coordinates by Damour, 
Schafer and Wex [6l|, [l^l and completed to 3PN order by Memmesheimer et al. [5^ . This 3PN parametrization gives 
the relative separation vector r = (r sin 0, r cos 0, 0) of the compact objects and the mean anomaly / as 

r = ar{l — Cr cosu), (2.1a) 

'^"^ = w + (/40 + /e^) sin 2t; + (.g40 + 560) sin 3?; 

+160 sin4t; + /i60 sin 5w, (2.1b) 



27r(t — to) . / \ / N 

/ = — = u~ etsm u+ [gu + get) [v - u) 



T 



+ iUt + kt) sin V + iat sin2w + h^t sin3w, (2.1c) 



where u is the eccentric anomaly, v = 2arctan{[(l + e^)/{l — e,/,)]^/^ tan(-u/2)} and T is the orbital period. The key 
element in the parametrization, that makes it useful for comparisons with numerical results, is that the auxiliary 
functions a^, e^, /40, /60, 540, 560, «60, ^60, n = 2tt/T, e*, 54*, 56t, fit, fet, iet, Kt and can be expressed 
exclusively in terms of the binding energy Eh, the total angular momentum L and the symmetric mass ratio 77 of the 
binary system. The complete expressions in terms of the dimensionless quantities E = E^/fJ- and h = L/{iiM) are 
listed in Eqs. (20) and Eqs. (25) of [s^ for ADM-type and harmonic coordinates, respectively. Here AI — Mi + M2 
and /i = MiMi/iMi + M2) are the total and reduced mass of the system, respectively. 
A comparison with the Newtonian accurate Keplerian parametrization 

r = a(l — ecosu), (2.2a) 

1 I \ 1/2 
1 + e \ u 
tan ■ 



4* ~ (t>() — 2arctan j ^^"^ 2' ' (2.2b) 

I = u — e sin u , (2.2c) 

illustrates that the concept of eccentricity is much more complex in general relativity and a single number, such as 
the Newtonian eccentricity e, no longer suffices to parametrize the shape of the orbit. Nevertheless, the similarity of 
the Newtonian and 3PN expressions suggest that the numbers et, and represent some measure of the deviation 
of the binary's orbit from quasi-circularity. This becomes particularly clear if we plot these quantities as functions of 
the orbital angular momentum L/M"^ for fixed binding energy E^/M and mass ratio 77. 

The result obtained for our sequence 1 models is shown in the left panel of Fig. [T] Several features of this plot 
are noteworthy. First, all eccentricity parameters diverge in the limit of a head-on collision. This is an artifact of 
the appearance of l/{—2Eh?) terms in the PN expressions for et, e^ and e^ in Eqs. (20) and (25) of Ref. |53]. We 
note that the limit L ^ Q also plays a special role in the Newtonian case. The usual distinction between the range 
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< e < 1, corresponding to bound elliptic orbits, and e > 1, corresponding to unbound parabolic or hyperbolic 
trajectories, no longer applies in the case of vanishing angular momenta. Since in Newtonian theory 

= l + f{A'h,M2)Ei,L^, (2.3) 

where /(Mi, M2) is a function of the masses, in the head-on limit we would formally have e = 1, irrespective of the 
sign of the binding energy. Indeed such trajectories only have one degree of freedom, the energy, and the concept of 
eccentricity no longer applies. In this sense, it is not surprising that the PN formalism fails to provide meaningful 
values for e*, and in the head-on limit. 

The second observation to be made is the steep gradient of all eccentricity parameters close to the circular limit 
of vanishing et, and e^. The strong sensitivity of these parameters to the orbital angular momentum L results 
in finite values of the three eccentricities (of about 0.1) even when using quasi-circular parameters, as obtained from 
Eq. (65) of Ref. [l6l |. Similarly, we observe that et, and do not vanish for the same values of L (see the inset 
in the left panel of Fig. [T|). Instead, the values of et and corresponding to the orbital angular momentum where 

vanishes are of the order of 0.1. A similar uncertainty results from comparing the PN values obtained in harmonic 
and Arnowitt-Deser-Misner- Transverse- Traceless (ADMTT) coordinates (cf. the results in the two gauges in the left 
panel). We thus take 0.1 as an approximate lower limit for these eccentricity parameters obtainable for such relatively 
large binding energies using the 3PN Keplerian parametrization. This is also approximately the value of et, and 
60 obtained for the quasi-circular configurations of Table HI 

A third noteworthy feature of the "quasi-Keplerian" PN parametrization (|2.ip is its breakdown for close, near- 
merger binary orbits. For example, if we tried to compute and for the "almost circular" parameters we use 
in this paper (as listed in Table |T|) they would turn out to be imaginary when, roughly speaking, P/M > 0.10 (or 
L/M'^ > 0.83). This is easy to understand by looking at the inset of the left panel of Fig. [T] There we see that these 
eccentricities have a zero crossing for values of L/M^ which are smaller than those specified in our quasi-circular 
simulations: in both ADMTT and harmonic coordinates, for the specified value of the binding energy goes to zero 
when L/M^ ~ 0.84, and goes to zero when L/Al'^ ~ 0.83. For L/M"^ larger than this "critical" value and 
become negative, so that Cr and are imaginary. In the case of sequences 2 and 3, we observe the same behaviour. 
This is just a sign that we should not trust the PN approximation for these highly relativistic configurations, so our 
eccentricity estimates should be taken with a grain of salt. 

An eccentricity plot using the binding energy of sequence 2 or 3 would look almost indistinguishable from the plot 
for sequence 1, as shown in the left panel of Fig. [1] so we decided not to display them here. Instead, in the right 
panel of Fig.[l]we show the eccentricities computed for a much smaller binding energy, E^/M — —0.001. This binding 
energy corresponds to a binary with much larger separation and smaller orbital velocity, that should be described 
with much higher accuracy by the quasi-Kcplerian PN parametrization. In fact, in the Newtonian limit the three 
eccentricities should agree with each other, reducing to the Newtonian definition at large separations. For example, 
to leading PN order et and are related by 

e^^et[l-{A-7i){2TTM/Tf\ (2.4) 

where T is the orbital period (see eg. [6^). The relation between the different eccentricities at higher PN orders can 
be found in Eq. (21) of Ref. [5|. 

From the right panel of Fig. [T] we see that the three eccentricity parameters do agree much better, as expected, when 
the binary members are far apart, and that differences resulting from the use of harmonic or ADMTT coordinates 
become negligible. We still see the breakdown of the formalism in the head-on limit. However, now all six curves are 
much closer to the expected Newtonian behavior, with vanishing eccentricity in the circular limit (where L approaches 
the maximum allowed value) and e « 1 for smaller angular momenta. Unfortunately, it is currently prohibitively costly 
from a computational point of view to start numerical simulations from such low binding energies. For this reason, 
in this paper we focus on the merger and ringdown signals resulting from eccentric binaries, rather than on detailed 
comparisons with PN predictions for the emission of GWs during the inspiral. 

B. Mora- Will parametrization 

An alternative estimate of the binary's initial eccentricity can be obtained using the PN diagnostic formalism 
developed by Mora and Will ([sO], henceforth MW). Instead of imposing a quasi-Keplerian parametrization with 
different eccentricities for i, r and 0, MW define a single eccentricity parameter cmw and a PN expansion parameter 
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C as follows: 



BMW 



c 



/^im^ + VMn^Y' ^^^^ 



Here flp and Qa are the orbital angular frequencies when the binary passes through a local maximum (pericenter) 
and through the next local minimum (apocenter), respectively. 

The MW eccentricity parameter gmw has several advantages: it is defined in terms of observable quantities, it 
reduces to its Newtonian counterpart for small orbital frequencies, and it is gauge invariant through first PN order. 
The binding energy and angular momentum of the system can be expressed as functions of gmw Sind The relevant 
equations for black hole binaries are given by Eqs. (2) and (3) in [63|. In Ref. these same PN equations have been 
used to study truly eccentric black hole binary initial data, and to point out some interesting features of the resulting 
effective potentials. 

The PN expansion parameter is related to the frequencies at periastron and apoastron by 

^ (l-eMW^/^ (l + eMw)4/3- ^ ■ ) 

To estimate the eccentricity, we can assume that the binary's orbit is (say) at apoastron, so that Q ^ Qa- Then 
the binding energy and angular momentum can be expressed as functions of bmw a-nd MQ. We can equate these 
functions to the binding energies and angular momenta listed in Table ID 

E^'^ieuw^Mn) = Eb, (2.7a) 
L^'^ieuw.Mn) = L. (2.7b) 

If our assumption is correct, and the orbit is indeed at apoastron, this system will have a solution (bmW: Mfl) with 
eccentricity cmw > 0. It should be obvious from Eq. (12. 6p that a solution with eMW < would simply correspond to 
the binary being at periastron. 

In Ref. [6^ it was shown that, when comparing numerically computed quasiequilibrium binary black hole initial 
data against the MW PN diagnostic predictions for a given orbital frequency Mf2, energies are usually in better 
agreement than angular momenta. Unfortunately, our initial data sets specify the orbital angular momentum and 
the binding energy, but not the orbital frequency. Therefore, to estimate the binary's eccentricity we adopt two 
different procedures. The first procedure is based on numerically solving the system (j2.7p for the two unknowns 
(gmWj Mn), which are both considered as free parameters. This procedure is quite general. However it can involve a 
significant amount of systematic error, since the angular momentum can deviate by a large amount from numerical 
initial data even for small eccentricities and/or large separations, and we are trying to solve for both energy and 
angular momentum simultaneously. 

The second procedure is based on an assumption of quasi-circularity, and it consists of two steps. We first assume 
that BMW = and solve Eq. (|2.7ap for AIV,. For sequence 1, the orbital frequency obtained by imposing E^^{eMW — 
0,Mn) = Eb is Mn = 0.0513; for sequ ence 2, we find Mfl = 0.0431; for sequence 3, we find = 0.0212. Then 
we substitute this value of Mil in Eq. (j2.7bp and solve for the eccentricity. We call this solution CMW.circ, because 
it is obtained assuming that the energy function is very close to the circular prediction, and that deviations of the 
eccentricity from zero only affect the angular momentum. 

Table U lists the eccentricities et obtained using the definition by Memmesheimer et al. and the two eccentricities 
estimated from the MW diagnostic. Dashed entries mean that no solutions to the system (|2.7p exist for physically 
relevant values of the parameters. The MW eccentricity estimate eMW is consistent with the quasi-Keplerian parameter 
Ct: in fact, the two definitions are roughly in agreement, within factors of order unity. It is also encouraging that 
the minima of et and bmw correspond to the same simulations along each sequence: these variables seem to provide 
a reasonably accurate measure of deviations from circularity, even for binaries which are very close to merger. The 
second procedure (by construction) should become inconsistent for large eccentricities, but eMW,circ is much closer 
to the eccentricity estimates 0.01) obtained in [63|, and measured in longer quasi-circular inspiral simulations by 
different groups 1^2^, i36„ §7,, £8,] . 

In the following we will use (somewhat arbitrarily) the parameter et obtained in harmonic coordinates as a measure 
of deviations from circularity of the initial binary configuration. We should still bear in mind that this parameter 
deviates significantly from the eccentricity e of a Keplerian ellipse, in particular for small angular momenta. 
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III. NUMERICAL SIMULATIONS 

The simulations presented in this work have been performed with the Lean code , which is based on the Cactus 
computational toolkit (69j . The code employs the so-called moving-puncture method, i.e. it evolves initial data of 
puncture type [T^I using the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formulation of the Einstein equations 
and gauge conditions which allow the black holes to move across the computational domain. In contrast to the 
original version of the BSSN system, where the conformal factor ip is evolved in the form of the logarithmic variable 
(j> = \ntfj, we here evolve the variable x = e"*"^, as originally introduced in Ref. These equations are numerically 
app roximated using fourth-order spatial discretization for sequences 1 and 2, and sixth-order spatial discretization 
rrii for sequence 3. Integration in time is performed with the method of lines using the fourth order accurate Runge- 
Kutta (RK) scheme. Using the notation of Table 1 in Ref. JJJ, this corresponds to schemes labelled RKx4 and RKxe, 
respectively. We note, however, that neither version of the code is genuinely fourth or sixth order accurate, as second 



order ingredients are used at the refinement and outer boundaries (see |14| for details). 

The Lean code facilitates use of mesh- refinement via the Carpet package [t^, It^ . Specifically, the computational 
grid is represented by a nested set of Cartesian boxes with resolution successively increasing by factors of two. The 
innermost refinement levels are typically split into two components centered around either hole and following the 
black hole motion. Once the black hole separation decreases below a threshold value, the two components merge into 
one centered around the origin. A more detailed description of the Lean code is given in Ref. hj. 

GWs are extracted using the Newman-Penrose formalism. The Newman-Penrose scalar \['4 is calculated using the 
electromagnetic decomposition of the Weyl tensor, and decomposed into contributions of different multipoles using 
spherical harmonics -2Yim of spin- weight —2 according to 

*4 = XI -s^^^V'fm- (3.1) 
e,m 

For all simulations presented in this work we extract GWs at different extraction radii r^x- We use the results obtained 
at different radii to estimate the uncertainty arising from the use of a finite extraction radius. 

The calculation of apparent horizons is performed using Thornburg's AHFinderDirect [tI, [t^ and plays an 
important role in the construction of initial data sequences. 

We use initial data of puncture type as provided by Ansorg's TwoPuncture thorn JT^] in all simulations. For 
zero spins, an initial data set is uniquely determined by the bare mass parameters mi and m2 of the holes, their 
initial coordinate separation D and the Bowen-York [t^] parameters Pi and P2 for the black holes' linear momenta. 
Without loss of generality we always set P = Pi = — P2. The initial orbital angular momentum is then given by 
L = D X P. Because the black holes are initially located on the x-axis and orbit in the xy-plane, the initial angular 
momentum is given by its z-component, which we can write as L = = DP (where P = Py). 

We conclude this brief description of the numerical code with a summary of the variables used in the remainder 
of this work. We denote by Mi and M2 the initial black hole masses. The total mass is the sum of the individual 
masses, M = Mi + M2, and the reduced mass is = M1M2/M. We measure the mass ratio using either q = M1/M2 
or the "symmetric mass ratio" parameter -q = ji/M often used in PN studies. Finally, we denote the Arnowitt-Deser- 
Misner (ADM) mass by Madm, which is close to (but in general different from) the total black hole mass M . The 
normalization of dimensional quantities is performed as follows. Initial parameters are normalized using the total 
black hole mass M for compatibility with the PN relations. In contrast, we measure radiated energy and momenta 
(as well as time t and radii r) in units of the ADM mass, as is common in numerical studies. All normalizations will 
be clear from the use of M or A/adm in column or axis labels; the relation between the two normalizations can be 
worked out from Eq. p.2p below and from the binding energies of the two sequences we consider, which are given in 
the caption of Table [J 



A. Sequence of initial data 



The construction of a sequence of non-spinning, equal-mass binary initial configurations with constant binding 
energy is a two-step process. First, we determine a quasi-circular configuration using the 3PN accurate relation (65) 
of Ref. [iBl, which provides the initial linear momentum P/M for given black hole masses Mi and M2 and separation 
D jM. Without loss of generality we fix the scaling freedom of the numerical coordinates by setting Mi = M2 = 0.5. 
For fixed values of P and Z), there thus remains the task of finding bare mass parameters mi = m2 which result 
in black-hole masses of Mi = M2 = 0.5. This is done iteratively by using the Newton-Raphson method with initial 
guesses m^ — Mi. In the absence of spins, the black-hole masses are given by the irreducible masses, as calculated 
by the apparent horizon finder AHFinderDirect. For the quasi-circular model we calculate the binding energy 
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TABLE I: Sequences of models studied in this work. The irreducible mass of the individual black holes is Mi = M2 — 0.5 

in all simulations, and consequently /i = 1/4. The binding energy for sequence 1 (top 16 models in the table) is E^^^ /M — 

—0.014465, that for sequence 2 (following 8 models) is E^^^ /M = -0.013229 and that for sequence 3 (following 9 models) is 
(3) 

E^ = —0.008861. The bottom line lists parameters for a model which serves for estimating the uncertainties of the long 
simulations of sequence 3 (see Sec. HIT B|l and has merely been included for completeness. AtcAu is the time of formation of a 
common apparent horizon as measured from the radiation peak, shifted by the extraction radius: Ate ah = ^cah + '"ex — tpeak- 
For sequence 3, we have not calculated apparent horizons and we estimate this quantity to be AfcAH = — 16 A4"adm. A'waves 



and Npunc are the number of orbits obtained from the phase of the {£ = 2 



2) multipole and the puncture's orbital motion. 



respectively. The last four columns list various PN estimates of the eccentricity (see Sec. |lT]for a discussion). 
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according to 

= AfADM - M. (3.2) 

The second step consists in the construction of additional configurations with the same binding energy but different 
hnear momentum parameter P and, thus, different orbital angular momentum L = DP. For this purpose we fix P 
and demand, as before. Mi = M2 ~ 0.5. We then iteratively solve for the numerical parameters mi — m2 and D 
which yield the required black hole masses and binding energy. 

Three sequences obtained this way with E\,/M = -0.014465, -0.013229 and -0.008861 are listed in Table HI In 
the following we will refer to them as "sequence 1" , "sequence 2" and "sequence 3" , respectively. The corresponding 
quasi-circular configurations are those with linear momentum P/M — 0.1383, 0.1247 and 0.0850. 

It is worth mentioning, in this context, that our choice of fixing the binding energy is not unique. Alternative choices 
in constructing sequences include keeping constant the coordinate separations of the punctures, or the proper horizon- 
to-horizon distance. We have opted against these two options because they would result in very short merger times 
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for small angular momenta. At constant binding energy, instead, smaller values of the orbital angular momentum 
imply larger separations of the holes and, thus, a delay in the formation of the common apparent horizon. 

B. Accuracy of the simulations 

In order to estimate the uncertainties associated with the numerically calculated quantities of sequences 1 and 2, we 
have performed a convergence analysis for the quasi-circular configuration with P/M — 0.1247 of sequence 2, and for 
the eccentric configuration with P/M — 0.08 of sequence 1. We have evolved these systems using the following grid 
setup: {(192, 96, 56, 24, 12) x (3, 1.5, 0.75) , hi} . By this notation we mean that there is a total of 8 refinement levels. 
The 5 outer levels are centered on the origin and extend out io x = y = z = ±192, 96, 56, 24 and 12 respectively. The 
3 inner levels have 2 components each, centered on either hole with radius 3, 1.5 and 0.75. Finally, the grid spacing 
is hi on the finest level (where hi — 1/48, /12 — 1/44 and /13 = 1/40) and increases by a factor of 2 consecutively on 
each level. 

Fourth-order convergence is shown in Fig. [5] for the (£ = 2, m = 2) multipole of the Newman-Penrose scalar 5*4, 
the total radiated energy E and the radiated angular momentum in the z-direction Jrad- 

Using the fourth order convergence, we apply a Richardson extrapolation to the total radiated energy and obtain 
E/Mp^Y^M = 0.03686, 0.03668 and 0.03656 respectively at extraction radii Tox/Madm = 50.7, 60.8 and 70.9. These 
values correspond very well to a 1/ rox fall-off of the uncertainty arising from the use of finite extraction radii. The total 
radiated energy extrapolated to — ^ 00 is E/Maum = 0.03583. For the medium resolution case (h = /i2 = 1/44) 
and using an extraction radius Tcx/Mabm = 70.9, this analysis predicts an uncertainty ~ 2% due to the discretization 
and ~ 2.5% due to the use of finite extraction radius. The uncertainties in the radiated angular momentum Jrad are 
~ 2% for both error sources. The convergence study of the eccentric simulation with P/M = 0.08 of sequence 1 yields 
similar error estimates. We estimate the resulting total error from quadratic error propagation to be about 3%. In 
fact, this estimate is likely to be very conservative because the two error sources have opposite signs: finite resolution 
tends to underestimate radiated energy and momenta, while finite extraction radius usually leads to an overestimate. 

Performing a convergence analysis of simulations lasting as long as those of sequence 3 requires vast computational 
resources. In order to reduce the cost, we view these simulations as part of a wider parameter study, to be presented 
elsewhere, which also involves unequal-mass binaries. In order to assess the accuracy of those long simulations, 
we have focussed on a quasi-circular configuration with q — 2 which is listed as the last entry in Table [D From 
experience, we consider unequal-mass binaries significantly more challenging numerically than systems with equal 
mass, and therefore feel justified in using the uncertainties resulting from this model as conservative error estimates 
of our sequence 3 runs. We have evolved this configuration using a grid {(384, 192, 128, 48, 24) x (6, 3, 1.5, 0.75), /i^} 
with hi = 1/44, 1/48 and 1/56. The bottom panels of Fig. [2] show the convergence of the phase and amplitude of ^'4 
(bottom right panel) and the radiated energy and angular momentum (bottom left panel). The difference between 
the higher resolution runs has been rescaled for second order convergence. We observe second order convergence for 
these runs except for the late stages near merger, when the convergence increases to third order. Similar glitches in 
the convergence near merger have been observed in [2^ . Using the same technique as above, we obtain uncertainties 
of about 3 % for the radiated energy and 6 % for the radiated angular momentum for the simulations using /12 = 1/48 
and rox = 80.7 Madm- 

Because most simulations have been performed at the medium resolution only, we cannot in general apply Richard- 
son extrapolation. Unless specified otherwise we therefore present results as obtained numerically using h = I /44 and 
Tox/Madm — 70.9 for sequences 1 and 2, as well a.s h — 1/48 and rox = 120 M for sequence 3, bearing in mind the 
uncertainties we have just mentioned. 

We conclude this discussion by mentioning that the Lean code has been demonstrated to accurately evolve black 
holes with large spins in Ref. [7^. Specifically, methods to calculate the spin from quasinormal ringing, balance 
arguments and apparent horizon data were shown to result in excellent agreement for Kerr parameters above 0.9. 

C. Numerical waveforms 

In Ref. [2I] we studied the multipolar energy distribution of unequal-mass black hole binaries in quasi-circular orbits. 
By projecting 2.5PN calculations of the inspiral gravitational waveforms onto spin-weighted spherical harmonics 
-2Yirm we concluded that odd-m multipoles of the radiation are suppressed for equal-mass binaries. Low-€ multipoles 
carry most of the radiation, and within a given /-multiplet, modes with £ = m are typically dominant for quasi-circular 
motion. This analytical prediction was shown to agree very well with numerical simulations, and it has recently been 
confirmed by more accurate PN calculations 't^. Since in this paper we study equal-mass runs, we expect even-m. 
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FIG. 2: Upper panels: convergence analysis of the P/M — 0.1247, quasi-circular model of sequence 2, using resolutions 
h = 1/48, 1/44 and 1/40. The left panel shows the differences in the {£ = 2, m = 2) multipole of Madm'"'i/'22, the right 
panel the total radiated energy and z-component of the angular momentum. In both cases, the differences between the higher 
resolution simulations have been rescaled by Q4 for the expected fourth order convergence. Middle panels: same, but for the 
P/M = 0.08, eccentric model of sequence 1. Lower panels: same, but for the unequal mass model listed at the bottom of Table 
[T] For clarity, we present the convergence of 'I'4 using phase and amplitude instead of the real part. For this configuration, 
we observe second order convergence. The corresponding convergence factors for the resolutions used here, are Qi = 1.58 and 
Q2 = 0.72. 



low-£ multipoles to be dominant, at least when the eccentricity is small enough. This expectation is again confirmed 
by our numerical time-evolutions. 

To be more quantitative, in Fig. [3] we show the modulus of the real part of the {£ = 2, m = 2), {£ = 2, m = 0) 
and (£ = 4, m = 4) components of the Weyl scalar for some representative runs. From these plots it is clear that the 
{£ = 2, m — 0) component contributes significantly for P/M < 0.05. That the {£ = 2, ni = 0) component should be 
significantly excited in the head-on limit P/M is known from previous numerical simulations: in fact, the m = 
mode would by far be dominant if we had chosen the collision to happen along the z-axis (see eg. [13, [8l|). Notice 
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FIG. 3: Modulus of the real part of the {£ — 2, m = 2), {£ = 2, m — 0) and {1 = 4:, m = 4) components of the waveforms. 
We show waveforms from four representative simulations: sequence 2 runs with P/M = (head-on limit), P/M = 0.02, 
P/M = 0.04 and P/M = 0.1247 (quasi-circular case). 

that in this paper GWs have always been extracted assuming that the z-axis is perpendicular to the orbital plane. In 
the case of head-on collisions this is contrary to most previous studies, where to take full advantage of the symmetry 
of the problem the axis of reference for the angular coordinates is identified with the axis of collision. In consequence, 
we find the radiated energy of head-on collisions to be quadrupole dominated, but to contain m — ±2 and m = 
contributions of comparable magnitude rather than almost exclusively an £ — 2, m — contribution as is the case 
for alignment of the two axes. Our choice is entirely motivated by using identical angular coordinates throughout the 
sequence of models. A detailed analysis of the transformation properties of multipolar components of the radiation 
under rotations, translations and boosts can be found in Ref. [82]. 

As P/M increases and we approach quasi-circular motion, the low-amplitude portion of the {£ = 2, m = 0) rnode 
decreases in amplitude, and it is significantly contaminated by noise. As expected from our previous study [2J|, in 
the same limit the amplitude of the {£ = A, m = 4) mode grows. Unfortunately, Fig. [3] illustrates that even for 
{£ = A, m = 4), which is the largest of the subdominant radiation modes, the ringdown signal is strongly distorted by 
either non-linear effects or boundary reflection noise. For this reason it is difficult to use higher multipoles to improve 
spin estimates from quasinormal mode (QNM) fittings, as proposed in [24]. This problem will be discussed in more 
detail in Sec. IVII below. 

We can obtain an estimate of the number of orbits in our simulations from the puncture trajectories, calculated 
according to dx'^/dt = —f3'^- For illustration, in Fig. [5] we plot the trajectories of four models of sequence 2 with 
P/M — 0.02, 0.06, 0.10 and 0.1247. The figure demonstrates the inspiralling nature of the simulations with large 
initial momentum, whereas those with small momentum rather represent plunging configurations. 

We define a phase 0punc of these trajectories by expressing the puncture's position in spherical polar coordinates 

(a;,y) = (rpuncCOS^punc, ^-punc sin 0punc)- (3.3) 

Then we consider the phase difference A^punc = </'punc(icah) — 0punc(^o)i where icah is the time of formation of a 
common apparent horizon and we choose to = SOA/adm to cut off the initial data burst (this is consistent with 
the derivation of the frequency from the GW signal, discussed below). The number of orbits then follows from 
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FIG. 4; Trajectories of the models with P/M = 0.02, 0.06, 0.10 and 0.1247 of sequence 2. The trajectory of one hole only is 
shown in each case. The positions of the respective second holes follow from symmetry across the origin. The + denote the 
locations of the individual holes at the time of common apparent horizon formation. 



^punc = A0punc/(27r). 

The gravitational wave signal \E'4 also serves to estimate the number of orbits completed by the binary prior to 
merger. For this purpose we focus on the (£ = 2, m = 2) multipole and decompose the mode coefficient into a 
time-dependent amplitude A{t) and phase (/)(t) according to 

V'22W =v4(<)e''>(*). (3.4) 

We next calculate the phase difference 

= (j){tcah + ''ox) - 0(^0 + ?'cx) , (3.5) 

where takes (approximately) into account the time it takes for the waves to propagate to the extraction radius. 
The number of orbits completed by the binary is then estimated as iV^avos = A(/)/(47r), where the additional factor of 
2 compensates for the difference between the orbital frequency and that of a multipole with to = 2. Both estimates 
of the number of orbits, together with the time of formation of the common apparent horizon, are given in Table [H 
Small differences in these numbers are due to the fact that the approximate relation cjwavcs — rnuo-p-anc breaks down 
near the merger of the holes. 

This is illustrated in Fig.[5l There we plot the frequency 0^22 and the phase (j)22 obtained from -022, comparing with 
the frequency mwpunc and phase (/)punc obtained from the puncture trajectory (see [2^ for details). Estimates from 
the punctures' motion are physically irrelevant after formation of the apparent horizon. The upper panels refer to 
different models from sequence 2: a nearly plunging motion with P/M — 0.06 (left) and a quasi-circular orbit with 
P/M = 0.1247 (right). It is clear from the figure that deviations between cjwavos and TOWpunc grow significantly near 
merger. The same holds for the two models of sequence 3 shown in the bottom panels. 

From Fig. [H and from the data for TVpunc and iVwavcs listed in Table HI we deduce that the simulations with L < 
Lcrit — 0.08M^ complete significantly less than one cycle, so they are effectively plunging trajectories. Consistently 
with this interpretation, we have see in Sec. HIl below that PN estimates of the eccentricity of the orbit fail for these 
plunging configurations. 



D. Polarization 



In Appendix D of [2J] we proposed to measure the polarization of a waveform using an "elliptical component of 
polarization" Pg. This quantity has the property that Pe = 1 for circular polarization, and Pe ~ for linear 
polarization. By looking at the dominant (£ — 2, to = 2) component of the radiation we also showed that, with the 
exception of the (unphysical) initial data burst and of the final part of the signal, which is dominated by noise, the 
polarization of a binary moving in a quasi-circular orbit is circular to a very good approximation. 

In Fig. [51 to be compared with Fig. 28 of Ref. [g^l, we show the real and imaginary parts of the dominant 
(£ = 2, TO = 2) multipolar component of the radiation emitted by equal-mass binaries with different values of P/M as 
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FIG. 5: Frequency and phase obtained from the £ = 2, m = 2 multipole of the gravitational radiation as well as the puncture 
trajectory for models P/AI = 0.06 and 0.1247 of sequence 2 (upper panels) and P/M = 0.06 and 0.0850 of sequence 3 (lower 
panels). The dotted vertical lines mark the formation of the apparent horizon. 



obtained from our sequence 2. In the bottom panel of each plot we compute the elliptical component of polarization. 
The results clearly illustrate that the polarization is linear in the head-on limit, where the imaginary component of 
the radiation vanishes, and that Pe ^ 1 as the orbit becomes circular. It is also worth noticing that the ringdown 
part of the signal is circularly polarized even when Pe is slightly less than one in the inspiral part (see eg. the plot 
for P/M = 0.02). 



IV. RADIATED ENERGY AND FINAL ANGULAR MOMENTUM 



In this section we study the multipolar energy distribution and the final angular momentum for the three sequences 
of simulations listed in Table |T1 

The multipolar distribution of radiated energy for sequence 1 is shown in the upper left panel of Fig. [71 We 
have excluded radiation due to the initial burst by ignoring contributions at simulation times i < 50 M + rex. The 
results demonstrate a relatively weak dependence of the radiated energy in each multipole on initial linear momenta 
P/M > 0.08, corresponding to an angular momentum of L/M"^ 0.8. 

For smaller linear (or angular) momenta, we observe that: (1) the total radiated energy decreases almost expo- 
nentially, (2) the relative contribution of multipoles with £ > 2 becomes weaker, and (3) the contribution of the 
(£ = 2, 171 = 0) mode increases to approximately the same level as the {£ = 2, m = ±2) modes. All of these features 
are quite insensitive to the inclusion of the spurious initial wave burst: only the (£ = 2, to = 0) mode is significantly 
contaminated by the initial radiation when P/M > 0.08. The same observations also apply to the models of se- 
quence 2, shown in the right panel of Fig. [71 Here the transition seems to occur at a linear momentum slightly below 
P/M = 0.08. From Table [J we see that this again corresponds to an angular momentum of L/M^ w 0.8. Similarly, 
the transition occurs just below P = 0.04 in the case of sequence 3, which again corresponds to an orbital angular 
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FIG. 6: i = m = 2 components of the low-L (or low-P) waveforms and their polarization for various models of sequence 2. 
For L = the imaginary part of the waveform is zero within the noise level (i.e., the cross component is zero for symmetry 
reasons). As L increases, the polarization becomes circular. Spikes in Pe at early times are due to the inital data burst, and 
spikes at late times are due to boundary reflection noise and low strength of the signal. 



momentum L/A/P « 0.8. A fit of sequence 1 runs by a polynomial in L/AP yields 

A similar behavior is found for the angular momentum of the final black hole, shown in Fig. [H] We have calculated 
the final spin from balance arguments: the final black hole mass Mgn is obtained by subtracting the total radiated 
energy from the initial ADM mass, and the final black hole angular momentum J^n is similarly given by the initial 
orbital angular momentum minus the momentum radiated in GWs. The figure shows the dimensionless Kerr parameter 
jfin — Jfin/M"^^. Again, small increases in eccentricity only lead to a mild increase in the final spin. The Kerr parameter 
has a local maximum, and then it decreases rapidly for LjlvP < 0.8. The results for the three sequences show 
remarkably good agreement below L/AP < 0.8 and merely differ at large angular momenta, as initial configurations 
with such large L only exist when choosing sufficiently large separations. By experimenting with sequence 1 runs, we 
found the following, reasonably accurate three-parameter polynomial fit: 

4^ = 0.0225 0.0381 f A) + 0.5589 f A) . (4.2) 



The large exponent in the final term is an artifact of the phenomenological nature of the fit and yields optimal 
agreement with the data. We found that the following alternative fitting functions, with two unknowns, also perform 
very well, yielding Er^d/MADM = 3.18 x 10"^ x e^-^"^^/^' and Jrad/M|j3M = 7-34 x IQ-^'e^-^''^/*^. 

An understanding of the apparent threshold in the orbital angular momentum L/M^ separating inspiraling and 
plunging orbits can be obtained by considering geodesic orbits in a Schwarzschild background. The key observation 
here is that particles of mass nip in closed orbits must satisfy the condition L/nip > 2\/3A4 . A rough extrapolation of 
this result to comparable-mass binaries would yield L/i-qM) > 2\/3M, or L/M^ > 0.866. There is no firm theoretical 
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P/M 



P/M 




FIG. 7: Multipolar energy distribution as a function of the initial momentum P for sequence 1 (upper left), sequence 2 (upper 
right) and sequence 3 (lower). Ee denotes the energy radiated in all multipoles with indices £ and m = —£, ...,£. We remove 
the initial data burst by ignoring all data with t < r^xt + 50M. 



justification for this extrapolation, and yet the point-particle threshold L/M"^ = 0.866 is remarkably close to the 
observed transition range of L/M"^ « 0.8. 

The agreement gets even better if, in the spirit of Ref. we use a slightly improved perturbative model, considering 
particle orbits around a Kerr black hole with spin given by the final Kerr parameter of our simulations (see SecfVland 
Appendix |X] for details). For eccentric inspirals, the minimum allowed angular momentum should be attained at the 
so-called separatrix corresponding to the maximal allowed eccentricity, Cp = 1, where is the eccentricity defined in 




FIG. 8: Final angular momentum as a function of the initial momentum P for sequence 1 (black solid, cross), 2 (red dashed, 
plus) and 3 (blue dotted curve, diamond). 
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Eq. (jA6p for point particle geodesies. The orbital angular momentum at the separatrix is also a function of the final 
black hole's spin, iscpOfim Cp), and it is defined below in Eq. (|5.ip . An explicit calculation yields striking agreement 
with numerical simulations: LscpOfin = 0.69, Cp = 1) = 0.778M^ when we consider the final spin j^n = 0.69 produced 
by a quasi-circular inspiral, and iscp(jfin = 0.724, Cp ~ 1) = 0.763M^ when the final hole has the largest spin observed 
in our simulations. 

Our interpretation of this similarity is that simulations with L/IvP < 0.8 can no longer be viewed as eccentric 
orbiting black-hole binaries, but rather represent plunging configurations or grazing collisions. The change in character 
of the energy distribution visible in Fig.[7]when L/M"^ w 0.8 thus demonstrates the transition from orbiting to plunging 
binary systems. 

Our simulations starting at different initial separation show some degree of universality in this transition. This 
can be understood in terms of the relatively low amount of angular momentum the binaries emit in the earlier 
stages of the inspiral (or plunge). In fact, consider some arbitrary configuration of sequence 3 with binding energy 

jrt\ 

/M = —0.008861. As the binary evolves, it emits gravitational radiation and the binding energy decreases. 
Of particular interest for our argument is the moment in time when the binding energy decreases to the value 

/OA 

El '/M = -0.013229 corresponding to sequence 2. In order to estimate this time, we consider the radiated energy 
Ai?rad(^) measured at r^x, integrated up to time t. The transition time t^2 from sequence 3 to sequence 2 is then 
approximated by the relation 

A^rad (^32 + r,^) = - E^^^ . (4.3) 

We next calculate the total amount of angular momentum A J].ad/-A^ADM(^32 +^ox) radiated away by the system during 
its transition from sequence 3 to sequence 2. The same procedure is applied for the transition from sequence 2 to 1. 
Results are shown in Table |TT1 

By comparison with the initial orbital angular momentum L/M^ listed in the third column of Table U we see that 
all binary models starting with orbital angular momentum near the critical value L/M"^ ~ 0.8 radiate only a few 
per cent of their angular momentum until they reach a more strongly bound state. This explains the approximate 
universality of models belonging to different sequences, but having the same initial orbital angular momentum. 

V. PERTURBATIVE ESTIMATES OF THE FINAL ANGULAR MOMENTUM OF ECCENTRIC 

BINARY BLACK HOLE COALESCENCE 

A simple, surprisingly successful model to compute the spin of the final black hole for binaries in quasicircular 
orbits has recently been proposed by Buonanno, Kidder and Lehner Here we discuss extensions of that model to 
eccentric binaries. The model introduced in Ref. 2] is based on three main assumptions: (i) the gravitational energy 
radiated after the formation of a common apparent horizon is a second order quantity, and has a small effect on the 
spin of the final hole; (ii) the magnitude of the individual spins remains constant during the inspiral, and (iii) most 
of the angular momentum is radiated during the long inspiral stage until the system reaches the innermost stable 
circular orbit (ISCO). A crucial final ingredient of the model is to estimate the contribution of the orbital angular 
momentum to the final spin by associating it with the orbital angular momentum of a point particle orbiting a Kerr 
black hole with spin parameter jfin corresponding to the final hole. 

The generalization of this model to eccentric orbits is straightforward, if one identifies the ISCO with the separatrix, 
or innermost stable bound orbit [s^ls^l (see Appendix \K\ for details). Following in the general case of eccentric 



TABLE II: The transition times 1^2 (tii) it takes a sequence 3 (2) model to radiate an amount of gravitational wave energy 
corresponding to the difference in binding energy from sequence 2 (1). The amount of angular momentum AJrad/A^ADM for 
models near the critical orbital angular momentum Lcrit radiated up to this time is small compared with the initial orbital 
angular momentum L/M"^ listed in Table |l] 





seq3 seq2 




seq2 


seql 


L/M^ T32/MADM 


AJrad/MADM(*32) 


L/M'^ t2l/MADM 


AJrad/MADM(i2l) 


1.020 


1646 


0.145 


0.873 


91 


0.020 


1.021 


1653 


0.145 


0.873 


98 


0.022 


1.011 


1450 


0.130 


0.868 


106 


0.023 


0.984 


926 


0.112 


0.839 


106 


0.021 


0.931 


444 


0.082 


0.765 


108 


0.017 


0.841 


203 


0.055 


0.612 


114 


0.014 


0.704 


201 


0.042 


0.350 


127 


0.010 
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orbits we write the (approximate) conservation equation 



. _ ^scpO'fin, ep) ^ MiQi ^ M2a2 .(.-.s 



where at is the Kerr parameter of each hole, the eccentricity Cp for point particle geodesies is defined in Eq. (jA6p 
and isep(ifin, Sp) is the orbital angular momentum at the separatrix, computed within a point-particle framework. 
A prescription for the computation of this quantity is given in Appendix |^ Note that Cp is defined in terms of 
coordinate positions of the particle, and therefore there is no direct physical relation between Cp and the eccentricity 
quantifiers discussed above. Assuming without loss of generality that Mi > AI2, we can rewrite Eq. (|5.ip as 

= hs^i^^ + ^(1 + 7134^)2 + |(1 - ./T^r , (5.2) 

where ji = ai/Mi and 77, as usual, denotes the symmetric mass ratio. 

We will focus on the special case where the initial spins have equal magnitude, ji — j2- In Fig- El (left panel) we 
show the predicted spin of the final hole as a function of the symmetric mass ratio 77 when both initial spins are 
aligned with the orbital angular momentum. The effect of eccentricity is always mild, and for all values of the initial 
spin magnitude, eccentricity tends to increase the spin of the final hole. 

Ref. [4I also predicted that, for circular orbits, the final spin of the remnant should increase as the mass ratio q ^ 1 
(i.e., as 77 — > 1/4) for ji < jcrit, while it should decrease as g — > 1 if ji > Jcrit- They estimated jcrit ~ 0.948. At 
the critical value, any merger will leave the final spin essentially unchanged, irrespective of the mass ratio q. This 
expectation is confirmed by our calculations (see the inset of the left panel of Fig. [S]). However, since a non-zero 
eccentricity always tends to increase the final spin, jcrit grows with eccentricity. 



TABLE III: Estimated critical value of the Kerr parameter jcrit (see text) for selected values of the eccentricity Cp, as defined 
perturbatively (see Appendix 

ep = 0.0 ep = 0.1 ep = 0.2 ep = 0.5 ep = 0.9 
0-948 0.950 0.953 0.972 0.998~ 



Table IIIII shows the critical value of the initial Kerr parameter j"'' — ji — j2 as a function of the eccentricity 
(as defined perturbatively). For large eccentricities ep ~ 1 the critical Kerr parameter is very close to the maximum 
possible value, j"'' '--^ 1: in other words, for large eccentricities the final spin should always increase as we approach 
the equal-mass limit. Therefore we conjecture that the maximum spin jgn ~ 0.724 that we found in our sequence of 
equal-mass merger simulations should be an absolute upper limit on the spin that can be achieved as the end-product 
of non-spinning black hole mergers: unequal-mass mergers should always produce smaller final spins. 




FIG. 9: Left panel: predicted angular momentum by point-particle extrapolation of results from perturbation theory. Solid, 
dashed and dash-dotted lines assume a point-particle eccentricity ep = 0, 0.5 and 0.9, respectively (see Appendix |X] for a 
definition of this eccentricity parameter, not to be confused with the Newtonian eccentricity). Numbers next to each set of 
lines denote different values of ji — j\ = j2, and the initial spin on each hole is assumed to be aligned with the orbital 
angular momentum. Right panel: we assume the initial spin on each hole to be antialigned with the orbital angular momentum 
{ji = ii = ji < 0). For various eccentricities (as indicated in the legend) we show the functional dependence between mass 
ratio and ji needed to produce a final non-spinning hole. 
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An interesting question explored in concerns spin-flip configurations. Suppose that initially both black holes 
have equal Kerr parameters ji = ji = j2, and spins antialigned with respect to the orbital angular momentum. 
What is the value of the symmetric mass ratio t] for which the inspiral produces a Schwarzschild black hole, as a 
function of ji? As argued in [3] these "critical" configurations could be particularly interesting, since mild variations 
of the parameters around the critical values may produce interesting orbital dynamics (eg. spin flips) and complex 
gravitational waveforms. 

The critical curve predicted by the model has been shown to yield good agreement with numerical simulations, 
especially for spins aligned with the orbital angular momentum [85, 86]. In particular, our study 86] obtains for a 
mass ratio q = A and initial dimensionless Kerr parameters of the individual holes ji = —0.75, —0.8 and —0.85 a 
final Kerr parameter jfi„ = 0.0533, 0.0237 and —0.0038, respectively, thus bracketing the formation of a Schwarzschild 
hole. For comparison, for zero eccentricity Eq. (|5.2p yields ji = —0.815083, while the fitting formula in 85] yields 
ji — —0.823462. Details of this study, together with a multipolar analysis of several spinning configurations, are 
given in Ref. ^86(]. As shown in the right panel of Fig. [9l the effect of eccentricity on the location of these spin-flip 
configurations is very mild. For a flxed mass-ratio, we predict that the magnitude of the (antialigned) spins required 
to produce a Schwarzschild black hole as the result of the merger should increase with eccentricity. 

Considering the remarkable success of point-particle extrapolations, some words of caution are required. First of all, 
when we extrapolate results to comparable masses we should not attach any special meaning to the actual values of the 
point-particle eccentricity, since we do not know how that notion translates into the different FN definitions discussed 
in Section ini At best, we can only expect that the FN eccentricities are related to the perturbative eccentricity, as 
implicitly defined in Eq. (jA6[) . by some monotonically increasing functional dependence. 

Furthermore, we should not forget that the model proposed in is only an approximation: for example, for a 
merger of equal-mass, non-spinning holes the model predicts a Kerr parameter jgn = 0.663 for the final hole, to be 
compared with jfjn — 0.691 from our numerical simulations (see also [24|). 

Despite these ambiguities and limitations, the point-particle model can still make remarkable predictions. For 
example, we can ask the following question: if we consider a non-spinning black hole binary, how much can we 
increase the Kerr parameter of the final hole by setting the binary in an eccentric orbit? We showed in Fig. [5] that, 
according to point-particle extrapolations, the final Kerr parameter should increase monotonically with eccentricity. 
Despite the unclear meaning of the eccentricity parameter, it is reasonable to assume that the maximum increase in the 
final Kerr parameter for non-spinning binaries should be Ajgjj*'^ = jfin(ep = 1) — jfin(ep = 0) ~ 0.750 — 0.663 ~ 0.087. 
As the motion turns into a plunge, jun should again decrease. From the results listed in Table U we can read out 
the "true", general relativistic prediction for the maximal spin increase induced by orbital eccentricity: ^j^^^^ = 
0.724 — 0.691 ~ 0.033. The level of agreement with simple extrapolations from perturbation theory is, once again, 
really surprising. 

Finally, it is interesting to ask if we can use the eccentricity-induced increase of the final Kerr parameter to violate 
the cosmic censorship conjecture, by producing a final object with jgn > 1. In principle this could be possible: even 
for zero eccentricity, if we consider maximally spinning black holes (ji = j2 = 1) with spins aligned to the orbital 
angular momentum, according to the point-particle ex;trapolation the final spin resulting from an equal-mass merger 
is jfin = 0.959, which is very close to the Kerr bound t3]- In fact, extrapolation of point-particle results suggests that 
cosmic censorship will not be violated. From the previous discussion it should be clear that the most "dangerous" 
situation corresponds to maximal eccentricities (cp = 1). In this case, the extrapolation of point-particle results can 
be performed analytically. The expression for general initial spins and general mass ratio is too cumbersome to display 
here. If we consider equal-mass binaries, we get 



Ji+j2 , 3 + V9-4(ji+j2) 

Jfin = + g . (5.3) 

From this equation we are led to conjecture that Kerr black holes produced as the result of a merger always satisfy 
the Kerr bound. The bound is only saturated when Cp = I and ji = j2 = 1- 

If instead we drop the equal-mass assumption, but we consider ji = ji = j2 (and again = 1), we get the following 
expression for the final spin: 



jfin = J^ + 2f] [{1 - j, - 77) + v/(l-277)(l-j,)+772j . (5.4) 

This suggests that in the ji 1 limit, j'fin — > 1 for any mass ratio. 

Early numerical simulations of spinning black hole binary coalescence support our conjecture: quasi-circular black 
hole binaries with spins aligned to the orbital angular momentum tend to radiate excess angular momentum by 
undergoing a longer inspiral, or "orbital hang-up," and they never seem to produce naked singularities 78, 85, 87, 88]. 
It will be interesting to see if this conclusion holds true, as we predict, when the black hole binary is both spinning 
and highly eccentric. 
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VI. THE RINGDOWN WAVEFORM 



It is well known that frequencies and quality factors of ringdown waveforms encode information about the properties 
of the Kerr black hole produced as a result of the merger (see ^| and references therein). Here we use the matrix pencil 
method to estimate time variations in the frequencies and quality factors of the dominant multipole (£ = 2, m = 2) 
as we reduce the orbital angular momentum of the orbit. This and other methods have been discussed at length in 
Refs. [UjIsO], and we will not repeat that discussion here. An introduction to Prony methods to estimate parameters 
of complex exponentials in noise, in the context of black hole perturbation theory, can be found in Ref. [o^. To 
remove boundary reflection noise from our simulations, we use the procedure described in [13] for quasi-circular 
inspiral simulations. 
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FIG. 10: Prony estimates for the quality factor Q (left) and dimensionless oscillation frequency Mujr (right) of the fundamental 
mode with £ — m — 2. Different linestyles (see the legend in the left panel) refer to different models in sequence 1. 



In Fig. \TU\ we show the time dependence of the estimated (dimensionless) QNM frequency Muju and of the quality 
factor Q for the dominant multipolar component (£ — m = 2) of the radiation. All quantities are plotted as functions 
of the starting time for the fit as measured from the time ipeak corresponding to the maximum in the modulus 
of the waveform's amplitude |'022(ipeak)|- We see the same features we observed in [l^, Is^- After a short transient 
for ^0 — ^poak ^ lOM, frequencies and quality factors become roughly constant, except at very late times, where the 
ringdown signal is very small and noise contaminates the estimates. Frequency estimates are usually better than 
quality factor estimates, and temporal variations in both quantities (which are probably induced by noise in the 
simulations) become more significant for the shorter, plunging configurations with low values of P/M. 



0.8 



; l=m=2 




































■- P/M=0.02 








P/M=0.03 








-■ P/M=0.04 








■- P/M=0.05 








- ■ Quasi-circular 













""o 10 20 30 40 50 



FIG. 11: Angular momentum estimate jun from Prony fits. Horizontal lines are estimates of the final Kerr parameter from 
energy balance arguments (see Table |l]). 



The quality factor is a dimensionless quantity, and as such it can only depend on the dimensionless Kerr parameter 
of the final hole: Q = (5(jfin)- The Kerr parameter can then be obtained by inverting this relation. We perform the 
inversion by interpolating QNM tables ^^9^ . Results are shown in Fig. [TT] There we discard all points for which we 
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get unphysical black hole parameters (eg. points for which the final estimated black hole mass Mgn > Madm)- 

Unfortunately, estimates of the final angular momentum become less reliable as i/Af^ — > 0. The reason is apparent 
by looking, for example, at the right panel of Fig. 5 in Ref. [s^. As the spin of the final hole decreases, and in 
particular for jfjn < 0.5, relatively large variations in the Kerr parameter produce mild variations in the £ — m — 2 
quality factor. Conversely, small oscillations in Q produce large variations in jfin- From Fig. [10] we see that noise- 
induced variations in Q actually increase for near head-on collisions. Therefore, QNM based estimates of the Kerr 
parameter are sensibly affected by noise when jgn ^ 0.5. The situation becomes even worse when we consider higher 
multipoles. However, it is clear from Fig. II II that QNM estimates are consistent (within the errors, which increase as 
P/M 0) with estimates obtained from angular momentum balance arguments. 



A. Energy-Maximized Orthogonal Projection 



As discussed thoroughly in [2j|, there is no unique definition of the ringdown starting time, and correspondingly 
there is no unique way to define the fraction of energy radiated into ringdown for a given waveform. A useful measure of 
the ringdown energy is provided by the "Energy Maximized Orthogonal Projection" (EMOP). In the EMOP was 
shown to yield reasonable estimates of the ringdown radiation from non-spinning, unequal-mass black hole binaries, 
and therefore we shall perform the same analysis on the present waveforms. 



TABLE IV: EMOP data for / = 2. Numbers separated by a comma correspond to the -I- and x polarizations, respectively. The 
fraction of the total energy in ringdown is strongly dependent on the eccentricity. We find that, independently of the run, the 
value of tfiMOP for a given polarization is generally at a fixed position relative to the maximum of the waveform's amplitude 
ipeak- We measure this relative difference by AfEMOP = ipcak — teMOP- Angular brackets denote an average over the two 
polarizations. 
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0.14 


0.42,0.41 


230.5 


9.1,6.1 
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0.0164 


0.1383 


0.40,0.41 
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4.9,8.4 
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0.0149 


0.13 


0.43,0.40 
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8.1,4.6 
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0.42,0.43 


222.2 
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0.09 
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0.07 
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0.66,0.75 
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4.3,8.3 


6.3 
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0.02 
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0.000475 



Our results are summarized in Table HVl for sequence 1. Results for sequence 2 and 3 are very similar, and we do 
not show them here. Ref. [l^l found that ringdown accounts for about 42% of the total energy emitted in the merger 
of non-spinning black hole binaries in quasi-circular orbits, i.e., -E'EMOp/^'rad — 42%. In the second column of Table 
IIVI we show the corresponding results for the present simulations. As a general trend, higher eccentricity means that 
the black holes spend less time orbiting each other before merger. Therefore one expects that the relative ringdown 
content, as measured by -EsMOpZ-^'rad, should increase with eccentricity, since less energy is radiated during the pre- 
merger phase. This trend is clearly visible in Table HVl On the other hand, the ringdown starting time (tEMOp)j 
measured from tpcak, is insensitive to the eccentricity of the run. This is consistent with the idea that ringdown is 
associated with the formation of a deformed common apparent horizon, and it should not depend on the details of 
the process leading to the formation of the horizon. Finally, in the last entry of Table IIVI we show the fraction of 
the total ADM mass radiated in ringdown waves. This number is important for estimates of ringdown detectability 
with both Earth-based GW detectors and LISA [s^, [9l| . A polynomial fit as a function of the dimensionless orbital 
angular momentum yields 



Mabm 



0.11 - 1.55 



L 



4.22 



L 



(6.1) 
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where (as in Table ITV| angular brackets denote an average over the two polarizations. 

VII. CONCLUSIONS AND OUTLOOK 

We have presented a study of the gravitational waveforms produced by sequences of equal-mass, non-spinning black 
hole binaries. For each sequence, the binding energy of the system is kept constant and the orbital angular momentum 
is progressively reduced to zero, producing orbits of increasing eccentricity and eventually a head-on collision. We 
find that the motion transitions from inspiral to plunge when the orbital angular momentum L = Lent — O-SM. 
For L < Lcrit the binary always completes less than 1 orbit, and PN estimates of the orbital eccentricity are no 
longer meaningful. As the initial momentum of the holes P/M the polarization quickly becomes linear, rather 
than circular, and the radiated energy drops (roughly) exponentially. For equal-mass, non-spinning binaries, orbits 
with L ~ Lcrit produce the largest dimensionless Kerr parameter for the remnant, jfin = J/AP ~ 0.724 ± 0.13 (to 
be compared with the Kerr parameter jfln — 0.69 resulting from quasi-circular inspirals). These results are quite 
insensitive to the initial separation of the holes, and they can be understood using extrapolations from black hole 
perturbation theory. Larger separations, as used in sequence 3, will be necessary to perform accurate comparisons 
with PN predictions for the evolution of eccentric binaries (see [12]). Such an analysis is beyond the scope of this 
work, however, and a corresponding analysis of the sequence 3 data will be presented elsewhere. 

For equal masses, we found that eccentric binary mergers with L ~ Lcmt maximize the Kerr parameter of the 
final black hole. Using arguments based on point-particle extrapolations, we proposed two conjectures: (1) jun — 
0.724±0.13 should be close to the largest Kerr parameter that can be produced by any non-spinning black hole binary 
merger, independently of the binary's mass ratio; and (2) even if we consider maximally spinning holes with spins 
aligned with the orbital angular momentum, orbital eccentricity should not lead to violations of the cosmic censorship 
conjecture. It will be interesting to check these conjectures using numerical simulations of eccentric binaries with 
unequal masses and non-zero spins. 

Acknowledgments 

We thank Clifford Will, Achamveedu Gopakumar, Christian Konigsdorffer and Gerhard Schafer for discussions 
about eccentricity in the PN formalism, and Luciano RezzoUa and Michele Vallisneri for comments on the manuscript. 
This work was supported in part by DFG grant SFB/Transregio 7 "Gravitational Wave Astronomy". E.B.'s research 
was supported by an appointment to the NASA Postdoctoral Program at the Jet Propulsion Laboratory, California 
Institute of Technology, administered by Oak Ridge Associated Universities through a contract with NASA; by the 
National Science Foundation, under grant number PHY 03-53180; and by NASA, under grant number NNG06GI60 to 
Washington University. V.C.'s work was partially funded by Fundagao para a Ciencia e Tecnologia (FCT) - Portugal 
through projects PTDC/FIS/64175/2006 and POCI/FP/81915/2007. We thank the DEISA Consortium (co-funded 
by the EU, FP6 project 508830), for support within the DEISA Extreme Computing Initiative (www.deisa.org). 
Computations were performed at LRZ Munich and the Doppler and Kepler clusters at the Institute of Theoretical 
Physics of the University of Jena. We are grateful to the Center for Computational Physics (CFC) in Coimbra for 
granting us access to the Milipeia cluster. 

APPENDIX A: ORBITS OF POINT PARTICLES IN THE KERR GEOMETRY AND THE ANGULAR 

MOMENTUM AT THE SEPARATRIX 

Here we consider a test body of mass moving in a Kerr spacetime, and we neglect radiation reaction (therefore 
we focus on purely geodesic motion). For simplicity, we only consider equatorial orbits. The equations of motion are 
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given in Boyer-Lindquist coordinates by 

,2 J = ±(v;)i/2 ^ (Ai) 

^ = V0 = -{MjunE - L) H — , (A2) 

r^^ = K ^ -MjUMj^r^E -L)+ (^' + M'jL)T ^^3^ 
or A 

0(r) = 7r/2 , (A4) 

where T = E{r^ + APjlJ - LMjen, Vr = - A(r^ + {L- Mjf^^Ef), A = - 2Mr + A'Pjl^. The two constants 
of motion 

E = E/mp, L-L/nip, (A5) 

denote the orbit's specific energy and z-component of angular momentum respectively. We have prograde (retrograde) 
orbits according to whether L > (< 0) . Bound orbits require < i? < 1. A general bound equatorial orbit can 
equivalently be described either by the constants E and L, or by a semi-latus rectum p and an eccentricity Cp (with 
< Bp < 1). The semi-latus rectum measures the size of the orbit, and the eccentricity measures the degree of 
non-circularity. We define these parameters in terms of the two turning points of the orbit (r^ is the periastron and 
ra the apastron): 

rp = , _f , ra = ^ ■ (A6) 

± |— Cp -L Gp 

The specific angular momentum and energy can be computed by solving the simultaneous equations for E and L 

-,1/2 



E = 



(A7) 



9 _ -N{p, ep)TAl^^{p, Cp) 
" ^ 2Fl^) ' ^^^^ 

where the upper (lower) sign corresponds to prograde (retrograde) motion, and we have defined 

X = L-MjfinE. (A9) 

The explicit form of the functions F, N and A^: is 

Fip,ep) = ^[p^-2M{3 + el)p^ + M^'S + elfp-4M^jlJl-el)% (AfO) 

N{p,ep) = ^[-Mp^ + M\i3 + el)-jlJp-M'ji^{l + 3el)], {All) 

A,ip,ep) - N^^AFM^M^-p/Mf. (A12) 

In general there are three different radial turning points (i.e., roots of 14 = 0): the periastron at r = Vp, the 
apastron at r = ra and a third root r = rs, defining a forbidden region. The case with rp = corresponds to a 
marginally stable orbit: once at the periastron, the particle will enter into a circular orbit of radius rp. This location 
has been called the innermost stable bound orbit ^41], and it is the generalization of the ISCO for general eccentric 
orbits. At this stage the orbit has become unstable, so that a slight inwards "push" will drive the particle to plunge 
into the black hole. Therefore, stable bound orbits should satisfy r^ < rp. This translates to the inequality 

a;2(l-Kep)(3-ep) </ . (Af3) 

The boundary curve defined by 

p2 ^a;2(l-|-ep)(3-ep) (A14) 
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defines the separatrix of bound orbits. 

Equations (|A7p - (|A14p allow one to compute the energy and orbital angular momentum at the separatrix. Strictly 
speaking, these equations are only valid for point particles. The substitution 




(A15) 



should provide a reasonable extrapolation to the general finite mass ratio case. 
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